Filtered back projection (FBP) algorithm for computer tomography

ABSTRACT

Reconstructing images of objects spirally scanned with two-dimensional detectors with a novel algorithm. The image reconstruction process is proven to create an exact image of the object under the ideal circumstances. The algorithm has an FBP (Filtered Back Projection) structure and works very efficiently. The algorithm uses less computer power and combines the benefits of Exact Algorithms and Approximate algorithms.

This is a Divisional of application Ser. No. 10/143,160 filed May 10,2002 now U.S. Pat. No. 6,574,299.

This invention relates to computer tomography, and in particular toprocesses and systems for reconstructing three dimensional images fromthe data obtained by a spiral scan, and this invention claims thebenefit of priority to U.S. Provisional Application No. 60/312,827 filedAug. 16, 2001.

BACKGROUND AND PRIOR ART

Over the last thirty years, computer tomography (CT) has gone from imagereconstruction based on scanning in a slice-by-slice process to spiralscanning. From the 1970s to 1980s the slice-by-slice scanning was used.In this mode the incremental motions of the patient on the table throughthe gantry and the gantry rotations were performed one after another.Since the patient was stationary during the gantry rotations, thetrajectory of the x-ray source around the patient was circular.Pre-selected slices through the patient have been reconstructed usingthe data obtained by such circular scans. From the mid 1980s to presentday, spiral type scanning has become the preferred process for datacollection in CT. Under spiral scanning a table with the patientcontinuously moves through the gantry that is continuously rotatingabout the table. At first, spiral scanning has used one-dimensionaldetectors, which receive data in one dimension (a single row ofdetectors). Later, two-dimensional detectors, where multiple rows (twoor more rows) of detectors sit next to one another, have beenintroduced. In CT there have been significant problems for imagereconstruction especially for two-dimensional detectors. In what followsthe data provided by the two-dimensional detectors will be referred toas cone-beam (CB) data or CB projections.

For three-dimensional (also known as volumetric) image reconstructionfrom the data provided by a spiral scan with two-dimensional detectors,there are two known groups of algorithms: Exact algorithms andApproximate algorithms, that each have known problems. Under idealcircumstances, exact algorithms can provide a replication of an exactimage. Thus, one should expect that exact algorithms would produceimages of good quality even under non-ideal (that is, realistic)circumstances. However, exact algorithms can be known to take many hoursto provide an image reconstruction, and can take up great amounts ofcomputer power when being used. These algorithms can require keepingconsiderable amounts of cone beam projections in memory. Additionally,some exact algorithms can require large detector arrays to be operableand can have limits on the size of the patient being scanned.

Approximate algorithms possess a filtered back projection (FBP)structure, so they can produce an image very efficiently and using lesscomputing power than Exact algorithms. However, even under the idealcircumstances they produce an approximate image that may be similar tobut still different from the exact image. In particular, Approximatealgorithms can create artifacts, which are false features in an image.Under certain circumstances these artifacts could be quite severe.

To date, there are no known algorithms that can combine the beneficialattributes of Exact and Approximate algorithms into a single algorithmthat is capable of replicating an exact image under the idealcircumstances, uses small amounts of computer power, and reconstructsthe exact images in an efficient manner (i.e., using the FBP structure).Here and everywhere below by the phrase that the algorithm of theinvention reconstructs an exact image we will mean that in theory thealgorithm is capable of reconstructing an exact image. Since in reallife any data contains noise and other imperfections, no algorithm iscapable of reconstructing an exact image.

Image reconstruction has been proposed in many U.S. patents. See forexample, U.S. Pat. Nos.: 5,663,995 and 5,706,325 and 5,784,481 and6,014,419 to Hu; 5,881,123 and 5,926,521 and 6,130,930 and 6,233,303 toTam; 5,960,055 to Samaresekera et al.; 5,995,580 to Schaller; 6,009,142to Sauer; 6,072,851 to Sivers; 6,173,032 to Besson; 6,198,789 to Dafni;6,215,841 and 6,266,388 to Hsieh. However, none of the patents overcomeall of the deficiencies to image reconstruction referenced above.

SUMMARY OF THE INVENTION

A primary objective of the invention is to provide an improved processand system for reconstructing images of objects that have been scannedin a spiral fashion with two-dimensional detectors.

A secondary objective of the invention is to provide an improved processand system for reconstructing images of spirally scanned objects that isknown to theoretically be able to reconstruct an exact image and not anapproximate image.

A third objective of the invention is to provide an improved process andsystem for reconstructing images of spirally scanned objects thatcreates an exact image in an efficient manner using a filtered backprojection (FBP) structure.

A fourth objective of the invention is to provide an improved processand system for reconstructing images of spirally scanned objects thatcreates an exact image with minimal computer power.

A fifth objective of the invention is to provide an improved process andsystem for reconstructing images of spirally scanned objects thatcreates an exact image with an FBP structure.

A sixth objective of the invention is to provide an improved process andsystem for reconstructing images of spirally scanned objects with largerpitch, leading to faster scans than previous techniques.

A seventh objective of the invention is to provide an improved processand system for reconstructing images, of spirally scanned objects whichtake less time than current techniques, thereby allowing use in everydayclinical applications.

An eighth objective of the invention is to provide an improved processand system for reconstructing images of spirally scanned objects that isCB projection driven allowing for the algorithm to work simultaneouslywith the CB data acquisition.

A ninth objective of the invention is to provide an improved process andsystem for reconstructing images of spirally scanned objects that doesnot requiring storage for numerous CB projections in computer memory.

A tenth objective of the invention is to provide an improved process andsystem for reconstructing images of spirally scanned objects that allowsfor almost real time imaging to occur where images are displayed as soonas a slice measurement is completed.

A first preferred embodiment of the invention uses a six overall stepprocess for reconstructing the image of an object under a spiral scan.In a first step a current CB projection is measured. Next, a family oflines is identified on a detector according to a novel algorithm. Next,a computation of derivatives between neighboring projections occurs andis followed by a convolution of the derivatives with a filter alonglines from the selected family of line. Next, using the filtered data,the image is updated by performing back projection. Finally, thepreceding steps are repeated for each CB projection until an entireobject has been scanned. This embodiment works with keeping several(approximately 2-4) CB projections in memory at a time and uses onefamily of lines.

For the second embodiment, the novel algorithm allows for one CBprojection to be kept in memory at a time and one family of lines isused.

For the third embodiment, two families of lines can be used incombination with either one CB projection in memory or with several CBprojections in memory.

Further objects and advantages of this invention will be apparent fromthe following detailed description of a presently preferred embodiment,which is illustrated schematically in the accompanying drawings.

BRIEF DESCRIPTION OF THE FIGURES

FIG. 1 shows a typical arrangement of a patient on a table that moveswithin a rotating gantry having an x-ray tube source and a detectorarray, where cone beam projection data sets are received by the x-raydetector, and an image reconstruction process takes place in a computerwith a display for the reconstructed image.

FIG. 2 shows an overview of the basic process steps of the invention.

FIG. 3 shows mathematical notations of the spiral scan about the objectbeing scanned.

FIG. 4 illustrates a PI segment of an individual image reconstructionpoint.

FIG. 5 illustrates a stereographic projection from the current sourceposition on to the detector plane used in the algorithm for theinvention.

FIG. 6 illustrates various lines and curves, such as boundaries, on thedetector plane.

FIG. 7 illustrates a family of lines used in the algorithm of theinvention.

FIG. 8 is a four substep flow chart for identifying the set of lines,which corresponds to step 20 of FIG. 2.

FIG. 9 is a seven substep flow chart for preparation for filtering,which corresponds to step 30 of FIG. 2.

FIG. 10 is a seven substep flow chart for filtering, which correspondsto step 40 of FIG. 2.

FIG. 11 is an eight substep flow chart for back projection, whichcorresponds to step 50 of FIG. 2.

FIG. 12 illustrates the first family of lines used in Embodiment Threeof the invention.

FIG. 13 illustrates the second family of lines used in Embodiment Threeof the invention.

DESCRIPTION OF THE PREFERRED EMBODIMENTS

Before explaining the disclosed embodiments of the present invention indetail it is to be understood that the invention is not limited in itsapplication to the details of the particular arrangement shown since theinvention is capable of other embodiments. Also, the terminology usedherein is for the purpose of description and not of limitation.

First Embodiment

FIG. 1 shows a typical arrangement of a patient on a table that moveswithin a rotating gantry having an x-ray tube source and a detectorarray, where CB projections are received by the x-ray detector, and animage reconstruction process takes place in a computer 4 with a display6 for displaying the reconstructed image. For the subject invention, thedetector array is a two-dimensional detector array. For example, thearray can include two, three or more rows of plural detectors in eachrow. If three rows are used with each row having ten detectors, then oneCB projection set would be thirty individual x-ray detections.

FIG. 2 shows an overview of the basic process steps of the inventionthat occur during the image reconstruction process occurring in thecomputer 4 using a first embodiment.

The first embodiment works with keeping several (approximately 2-4) CBprojections in computer memory at a time and uses one family of lines.

In the first step 10, a current CB projection set is taken. The nextstep 20 identifies a set of lines on a virtual x-ray detector arrayaccording to the novel algorithm, which will be explained later ingreater detail. In the given description of the algorithm it is assumedthat the detector array is flat, so the selected line can be a straighttilted line across the array.

The next step 30 is the preparation for the filtering step, whichincludes computations of the necessary derivative of the CB projectiondata for the selected lines.

The next step 40 is the convolution of the computed derivative (theprocessed CB data) with a filter along lines from the selected family oflines. This step can also be described as shift-invariant filtering ofthe derivative of the CB projection data. In the next step 50, the imageof the object being computed is updated by performing back projection.

The invention will now be described in more detail by first describingthe main inversion formula followed by the novel algorithm.

The spiral path C of the x-ray source is described by the followingequations and depicted in FIG. 3, which shows mathematical notations ofthe spiral scan about the object being scanned:

y ₁(s)=R cos(s), y ₂(s)=R sin(s), y ₃(s)=s(h/2π),  (1)

Here

s is a real parameter;

h is pitch of the spiral;

R is distance from the x-ray source to the isocenter.

The object being scanned is located inside an imaginary cylinder U ofradius r, r<R (see FIG. 3). Let ψ be a smooth function with theproperties

ψ(0)=0; 0<ψ′(t)<1, tε[0,2π].  (2)

Even though it is not necessary, we will assume in addition that

ψ′(0)=0.5; ψ^((2k+1))(0)=0, k≧1.  (3)

Here and in what follows we assume that s₀, s₁, and s₂ are alwaysrelated by

s ₁=ψ(s ₂ −s ₀)+s ₀ if s ₀ ≦s ₂ <s ₀+2π,  (4)

s ₁=ψ(s ₀ −s ₂)+s ₂ if s ₀−2π<s ₂ <s ₀.  (5)

Conditions (2) and (3) can be easily satisfied. One can take, forexample, ψ(t)=t/2, and this leads to

s ₁=(s ₀ +s ₂)/2, s ₀−2π<s ₂ <s ₀+2π.  (6)

Denote $\begin{matrix}{{{u\left( {s_{0},s_{2}} \right)} = {\frac{\left( {{y\left( s_{1} \right)} - {y\left( s_{0} \right)}} \right) \times \left( {{y\left( s_{2} \right)} - {y\left( s_{0} \right)}} \right)}{\left| {\left( {{y\left( s_{1} \right)} - {y\left( s_{0} \right)}} \right) \times \left( {{y\left( s_{2} \right)} - {y\left( s_{0} \right)}} \right)} \right|}{sgn}\quad \left( {s_{2} - s_{0}} \right)}}{\left. {{{if}{\quad \quad}0} <} \middle| {s_{2} - s_{0}} \middle| {< {2\pi}} \right.,}} & (7) \\{{u\left( {s_{0},s_{2}} \right)} = {{\frac{{\overset{.}{y}\left( s_{0} \right)} \times {\overset{¨}{y}\left( s_{0} \right)}}{\left| {{\overset{.}{y}\left( s_{0} \right)} \times {\overset{¨}{y}\left( s_{0} \right)}} \right|}\quad {if}\quad s_{2}} = {s_{0}.}}} & (8)\end{matrix}$

Here

y(s₀),y(s₁),y(s₂) are three points on the spiral related according to(4), (5);

u(s₀,s₂) is a unit vector perpendicular to the plane containing thepoints

y(s₀),y(s₁),y(s₂);

{dot over (y)}(s):=dy/ds;

ÿ(s):=d²y/ds².

It is known in the literature that any point strictly inside the spiralbelongs to one and only one PI segment (see P. E. Danielsson et al.“Towards exact reconstruction for helical cone-beam scanning of longobjects. A new detector arrangement and a new completeness condition” in“Proc. 1997 Meeting on Fully 3D Image Reconstruction in Radiology andNuclear Medicine (Pittsburgh)”, eds. D. W. Townsend and P. E. Kinahan,yr. 1997, pp. 141-144, and M. Defrise, F. Noo, and H. Kudo “A solutionto the long-object problem in helical cone-beam tomography”, Physics inMedicine and Biology, volume 45, yr. 2000, pp. 623-643). A PI segment isa segment of line endpoints of which are located on the spiral andseparated by less than one pitch in the axial direction (see FIG. 4).Let s=s_(b)(x) and s=s₁(x) denote values of the parameter correspondingto the endpoints of the PI segment containing a reconstruction point x.We will call I_(PI)(x):=[s_(b)(x),s_(t)(x)] the PI parametric interval.The part of the spiral corresponding to I_(PI)(x) will be denotedC_(PI)(x) (see FIG. 4 which illustrates a PI segment of an individualimage reconstruction point). Next we fix a reconstruction point x insidethe spiral and s₀εI_(PI)(x). Find s₂εI_(PI)(x) such that the planethrough y(s₀), y(s₂), and y(s₁(s₀,s₂)) contains x. More precisely, wehave to solve for s₂ the following equation

(x−y(s ₀))·u(s ₀ , s ₂)=0, s ₂ εI _(PI)(x).  (9)

Such s₂ exists, is unique, and depends smoothly on s₀. Therefore, thisconstruction defines s₂:=s₂(s₀,x) and, consequently,u(s₀,x):=u(s₀,s₂(s₀,x)). Equation (9) can be solved by a variety ofmethods, all known under the name “root finders”. The mainreconstruction formula now is as follows: $\begin{matrix}{{{f(x)} = {{{{- \frac{1}{2\pi^{2}}}{\int_{I_{PI}{(x)}}{\frac{1}{\left| {x - {y(s)}} \right|}{\int_{0}^{2\pi}{\frac{\partial}{\partial q}{D_{f}\left( {{y(q)},{\Theta \left( {s,x,\gamma} \right)}} \right)}}}}}}}_{q = s}\frac{d\quad \gamma}{\sin \quad \gamma}{s}}},} & (10)\end{matrix}$

where

f is the function representing the distribution of the x-ray attenuationcoefficient inside the object being scanned,

e(s,x)=β(s,x)×u(s,x),

x is the cross-product of two vectors,

Θ(s,x,γ):=cos γβ(s,x)+sin γe(s,x),

D_(f) is the cone beam transform of f: $\begin{matrix}{{{D_{f}\left( {y,\Theta} \right)} = {\int_{0}^{\infty}{{f\left( {y + {\Theta \quad t}} \right)}{t}}}},} & (11)\end{matrix}$

${\beta \left( {s,x} \right)}:=\frac{x - {y(s)}}{\left| {x - {y(s)}} \right|}$

is the unit vector from the focal point y(s) pointing towards thereconstruction point x.

The proof that equation (10) gives a theoretically exact image ispresented in A. Katsevich “Improved exact FBP algorithm for Spiral CT”submitted for possible publication to the journal “Advances in AppliedMathematics” in November 2001.

Now we describe an efficient (that is, of the FBP type) implementationof inversion formula (10). It is clear from (9) that s₂(s,x) actuallydepends only on s and β(s,x). Therefore, we can write

u(s,β):=u(s,s ₂(s,β)),e(s,β):=β×u(s,β),βεS ², $\begin{matrix}{{{{{\Psi \left( {s,\beta} \right)}:={\int_{0}^{2\pi}{\frac{\partial}{\partial q}{D_{f}\left( {{y(q)},{{\cos \quad {\gamma\beta}} + {\sin \quad \gamma \quad {e\left( {s,\beta} \right)}}}} \right)}}}}}_{q = s}\frac{1}{\sin \quad \gamma}{\gamma}},} & (12) \\{{f(x)} = {{- \frac{1}{2\pi^{2}}}{\int_{I_{PI}{(x)}}{\frac{1}{\left| {x - {y(s)}} \right|}{\Psi \left( {s,{\beta \left( {s,x} \right)}} \right)}{{s}.}}}}} & (13)\end{matrix}$

Here S² is the unit sphere.

To better understand equations (12), (13), we illustrate variousimportant features on the detector array. Suppose first that the x-raysource is fixed at y(s) for some sεI_(PI)(x). Project stereographicallythe upper and lower turns of the spiral onto the detector plane as shownin FIG. 5 which illustrates a stereographic projection from the currentsource position on to the detector plane used in the algorithm for theinvention.

Since the detector array rotates together with the source, the detectorplane depends on s and is denoted DP(s). It is assumed that DP(s) isparallel to the axis of the spiral and is tangent to the cylinder y₁²+y₂ ²=R² (cf. equation (1)) at the point opposite to the source. Thus,the distance between y(s) and the detector plane is 2R. In order to makeFIG. 5 more readable, the detector plane is drawn slightly away from thespiral. Introduce coordinates in the detector plane as follows. Let thed₁-axis be perpendicular to the axis of the spiral, and thed₂-axis—parallel to it. This gives the following parametric curves:$\begin{matrix}{{{d_{1}(q)} = {2R\frac{\sin \left( {q - s} \right)}{1 - {\cos \left( {q - s} \right)}}}},\quad {{d_{2}(q)} = {\frac{h}{\pi}\frac{q - s}{1 - {\cos \left( {q - s} \right)}}}},} & (14) \\{{\Delta \leq {q - s} \leq {{2\pi} - {\Delta \quad o\quad r\quad \Delta} - {2\pi}} \leq {q - s} \leq {- \Delta}},} & \quad\end{matrix}$

where Δ is determined by the radius r of the imaginary cylinder U insidewhich the patient is located (see FIG. 3): Δ=2 cos⁻¹(r/R). The top andbottom curves are denoted Γ_(top) and Γ_(bot), respectively (see FIG. 6which illustrates various lines and curves, such as boundaries, on thedetector plane). The common asymptote of Γ_(top) and Γ_(bot) is denotedL₀. Let {circumflex over (x)} denote the projection of x. SincesεI_(PI)(x), {circumflex over (x)} is projected into the area betweenΓ_(top) and Γ_(bot) (see FIG. 6). Fix s₂ε[s−2π+Δ,s+2π−Δ],s₂≠s, and letΠ(s₂) denote the plane through y(s),y(s₂), and y(s₁(s,s₂)). If s₂=s,Π(s₂) is determined by continuity and coincides with the plane throughy(s) and parallel to {dot over (y)}(s),ÿ(s). The family of lines L(s₂)obtained by intersecting Π(s₂) with the detector plane is shown in FIG.7.

By construction, given any xεU with β(s,x) parallel to Π(s₂) and suchthat {circumflex over (x)} appears to the left (right) of the point ofwhere the line L(s₂) intersects Γ_(top) (Γ_(bot)) for the first time if{circumflex over (x)} is above (below) L₀, s₂ used here is precisely thesame as s₂ found by solving (9). The condition that we have formulatedregarding the location of {circumflex over (x)} relative to s₂ and L₀guarantees that s₂εI_(PI)(x). Since e(s,β)·β=0,|e(s,β)|=1, we can write:

β=(cos θ, sin θ);e(s,β)=(−sin θ,cos θ);β,e(s,β)εΠ(s ₂).  (15)

Therefore, $\begin{matrix}{{{{{\Psi \left( {s,\beta} \right)} = {\int_{0}^{2\pi}{\frac{\partial}{\partial q}{D_{f}\left( {{y(q)},\left( {{\cos \quad \left( {\theta + \gamma} \right)},{\sin \quad \left( {\theta + \gamma} \right)}} \right)} \right)}}}}}_{q = s}\frac{1}{\sin \quad \gamma}{\gamma}},{\beta \in {{\Pi \left( s_{2} \right)}.}}} & (16)\end{matrix}$

Equation (16) is of convolution type and one application of Fast FourierTransform (FFT) gives values of Ψ(s,β) for all βεΠ(s₂) at once.

Equations (13) and (16) would represent that the resulting algorithm isof the FBP type. This means that processing of every CB projectionconsists of two steps. First, shift-invariant and x-independentfiltering along a family of lines on the detector is performed. Second,the result is back-projected to update the image matrix. The mainproperty of the back-projection step is that for any point {circumflexover (x)} on the detector the value obtained by filtering at {circumflexover (x)} is used for all points x on the line segment connecting thecurrent source position y(s) with {circumflex over (x)}. Since ∂/∂q in(16) is a local operation, each CB projection is stored in memory assoon as it has been acquired for a short period of time for computingthis derivative at a few nearby points and is never used later.

Now we describe the algorithm in detail following the six steps 10-60shown in FIG. 2.

Step 10. Load the current CB(cone beam) projection into computer memory.Suppose that the mid point of the CB projections currently stored inmemory is y(s₀). The detector plane corresponding to the x-ray sourcelocated at y(s₀) is denoted DP(s₀).

Step 20. FIG. 8 is a four substep flow chart for identifying the set oflines, which corresponds to step 20 of FIG. 2. Referring to FIG. 8, theset of lines can be selected by the following substeps 21, 22, 23 and24.

Step 21. Choose a discrete set of values of the parameter s₂ inside theinterval [s₀−2π+Δ,s₀+2π−Δ].

Step 22. For each selected s₂ compute the vector u(s₀,s₂) according toequations (7), (8).

Step 23. For each u(s₀,s₂) computed in Step 22 find a line which isobtained by intersecting the plane through y(s₀) and perpendicular tothe said vector u(s₀,s₂) with the detector plane DP(s₀).

Step 24. The collection of lines constructed in Step 23 is the requiredset of lines (see FIG. 7 which illustrates a family of lines used in thealgorithm of the invention).

Step 30. Preparation for filtering

FIG. 9 is a seven substep flow chart for preparation for filtering,which corresponds to step 30 of FIG. 2, which will now be described.

Step 31. Fix a line L(s₂) from the said set of lines obtained in Step20.

Step 32. Parameterize points on the said line by polar angle γ in theplane through y(s₀) and L(s₂).

Step 33. Choose a discrete set of equidistant values γ_(j) that will beused later for discrete filtering in Step 40.

Step 34. For each γ_(j) find the unit vector β_(j) which points fromy(s₀) towards the point on L(s₂) that corresponds to γ_(j).

Step 35. Using the CB projection data D_(f)(y(q),Θ) for a few values ofq close to s₀ find numerically the derivative (∂/∂q)D_(f)(y(q),Θ)|_(q=s) ₀ for all Θ=β_(j).

Step 36. Store the computed values of the derivative in computer memory.

Step 37. Repeat Steps 31-36 for all lines L(s₂) identified in Step 20.This way we will create the processed CB data Ψ(s₀,β_(j)) correspondingto the x-ray source located at y(s₀).

Step 40. Filtering

FIG. 10 is a seven substep flow chart for filtering, which correspondsto step 40 of FIG. 2, which will now be described.

Step 41. Fix a line from the said family of lines identified in Step 20.

Step 42. Compute FFT of the values of the said processed CB datacomputed in Step 30 along the said line.

Step 43. Compute FFT of the filter 1/sin γ

Step 44. Multiply FTT of the filter 1/sin γ(the result of Steps 43) andFTT of the values of the said processed CB data (the result of Steps42).

Step 45. Take the inverse FFT of the result of Step 44.

Step 46. Store the result of Step 45 in computer memory.

Step 47. Repeat Steps 41-46 for all lines in the said family of lines.This will give the filtered CB data Φ(s₀,β_(j)).

By itself the filtering step is well known in the field and can beimplemented, for example, as shown and described in U.S. Pat. No.5,881,123 to Tam, which is incorporated by reference.

Step 50. Back-projection

FIG. 11 is an eight substep flow chart for back projection, whichcorresponds to step 50 of FIG. 2, which will now be described.

Step 51. Fix a reconstruction point x, which represents a point insidethe patient where it is required to reconstruct the image.

Step 52. If s₀ belongs to I_(PI)(x), then the said filtered CB dataaffects the image at x and one performs Steps 53-58. If s₀ is not insidethe interval I_(PI)(x), then the said filtered CB data is not used forimage reconstruction at x. In this case go back to Step 51 and chooseanother reconstruction point.

Step 53. Find the projection {circumflex over (x)} of x onto thedetector plane DP(s₀) and the unit vector) β(s₀,x), which points fromy(s₀) towards x.

Step 54. Using equation (9) identify the lines from the said family oflines and points on the said lines that are close to the said projection{circumflex over (x)}. This will give a few values of Φ(s₀,β_(j)) forβ_(j) close to β(s₀,x).

Step 55. With interpolation estimate the value of Φ(s₀,β(s₀,x)) from thesaid values of Φ(s₀,β_(j)) for β_(j) close to β(s₀,x).

Step 56. Compute the contribution from the said filtered CB data to theimage being reconstructed at the point x by dividing Φ(s₀,β(s₀,x)) by−2π²|x−y(s₀)|.

Step 57. Add the said contribution to the image being reconstructed atthe point x according to a pre-selected scheme (for example, theTrapezoidal scheme) for approximate evaluation of the integral inequation (15).

Step 58. Go to Step 51 and choose a different reconstruction point x.

Step 60. Go to Step 10 (FIG. 2) and load the next CB projection intocomputer memory. The image can be displayed at all reconstruction pointsx for which the image reconstruction process has been completed (thatis, all the subsequent CB projections are not needed for reconstructingthe image at those points). Discard from the computer memory all the CBprojections that are not needed for image reconstruction at points wherethe image reconstruction process has not completed. The algorithmconcludes when the scan is finished or the image reconstruction processhas completed at all the required points.

Embodiment Two

For the second embodiment, one CB (cone beam) projection can be kept inmemory at a time and, as before, only one family of lines on thedetector is used. Now we describe the second embodiment of the inventionin detail. Integrating by parts with respect to s in equation (10) weobtain an inversion formula in which all the derivatives are performedwith respect to the angular variables. $\begin{matrix}\begin{matrix}{{f(x)} = {{- \frac{1}{2\pi^{2}}}\left\{ {\left\lbrack {\frac{1}{\left| {x - {y(s)}} \right|}{\int_{0}^{2\pi}{{D_{f}\left( {{y(s)},{\Theta \left( {s,x,\gamma} \right)}} \right)}\frac{\gamma}{\sin \quad \gamma}}}} \right\rbrack |_{s = {s_{b}{(x)}}}^{s = {s_{t}{(x)}}} -} \right.}} \\{{{\int_{I_{PI}{(x)}}{\left( {\frac{\partial}{\partial s}\frac{1}{\left| {x - {y(s)}} \right|}} \right){\int_{0}^{2\pi}{{D_{f}\left( {{y(s)},{\Theta \left( {s,x,\gamma} \right)}} \right)}\frac{\gamma}{\sin \quad \gamma}{s}}}}} -}} \\{{\int_{I_{PI}{(x)}}\frac{{\beta_{s}\left( {s,x} \right)} \cdot {u\left( {s,x} \right)}}{\left| {x - {y(s)}} \right|}}} \\{{{\int_{0}^{2\pi}{\left( {\nabla_{u{({s,x})}}D_{f}} \right)\left( {{y(s)},{\Theta \left( {s,x,\gamma} \right)}} \right)\quad {\cot (\gamma)}{\gamma}{s}}} -}} \\{{{\int_{I_{PI}{(x)}}{\frac{{e_{s}\left( {s,x} \right)} \cdot {u\left( {s,x} \right)}}{\left| {x - {y(s)}} \right|}{\int_{0}^{2\pi}{\left( {\nabla_{u{({s,x})}}D_{f}} \right)\left( {{y(s)},{\Theta \left( {s,x,\gamma} \right)}} \right){\gamma}{s}}}}} -}} \\{{\int_{I_{PI}{(x)}}\frac{{\beta_{s}\left( {s,x} \right)} \cdot {e\left( {s,x} \right)}}{\left| {x - {y(s)}} \right|}}} \\{\left. {\int_{0}^{2\pi}{\left( {\frac{\partial}{\partial\gamma}{D_{f}\left( {{y(s)},{\Theta \left( {s,x,\gamma} \right)}} \right)}} \right)\frac{\gamma}{\sin \quad \gamma}{s}}} \right\}.}\end{matrix} & (17)\end{matrix}$

Here

β_(s)=∂β/∂s;

e_(s)=∂e/∂s, and

∇_(u)D_(f) denotes the derivative of D_(f) with respect to the angularvariables along the direction of the vector u: $\begin{matrix}{{{{\left( {\nabla_{u}D_{f}} \right)\left( {{y(s)},\Theta} \right)} = {\frac{\partial}{\partial t}{D_{f}\left( {{y(s)},{{\sqrt{1 - t^{2}}\Theta} + {t\quad u}}} \right)}}}}_{t = 0},{\Theta \in {u^{\bot}.}}} & (18)\end{matrix}$

Comparing equation (10) and equation (17) we see that equation (17)admits absolutely analogous filtered back-projection implementation.Moreover, since no derivative with respect to the parameter along thespiral is present, there is never a need to keep more than one CBprojection in computer memory at a time. Now we describe the algorithmin detail.

Step 10. Load the current CB projection into computer memory and discardthe CB projection that was in computer memory before. Suppose that theCB projection just loaded into computer memory corresponds to the x-raysource located at y(s).

Step 20 is the same as Step 20 in Embodiment One with s₀ replaced by s.

Step 30. Preparation for filtering

Steps 31-34 are the same as in Embodiment One with s₀ replaced by s.

Step 35. Using the CB projection data D_(f)(y(s),Θ) findD_(f)(y(s),β_(j)) and the derivatives(∇_(u(s,x))D_(f))(y(s),Θ(s,x,γ))|_(Θ=β) _(j) and${{\frac{\partial}{\partial\gamma}{D_{f}\left( {{y(s)},{\Theta \left( {s,x,\gamma} \right)}} \right)}}}_{\Theta = \beta_{j}}$

 for all β_(j). This can be done using interpolation and finitedifferences.

Step 36. Store the values computed in Step 35 in computer memory.

Step 37. Repeat Steps 31-36 for all lines L(s₂) identified in Step 20.This way we will create the processed CB data D_(f)(y(s),β_(j)),(∇_(u(s,x))D_(f))(y(s),β_(j)), and$\frac{\partial\quad}{\partial\gamma}{{D_{f}\left( {{y(s)},\beta_{j}} \right)}.}$

Step 40. Filtering

Step 41. Fix a line from the said family of lines identified in Step 20.

Step 42. Using FFT convolve the said processed CB data computed in Step30 with filters 1/sin γ and cot γ along the said line according toequation (17). This will give the following three kinds of the filteredCB data (see also equations (12), (15), and (16)): $\begin{matrix}{{{\Psi_{1}\left( {s,\beta} \right)} = {\int_{0}^{2\quad \pi}{{D_{f}\left( {{y(s)},{{\cos \quad \gamma \quad \beta} + {\sin \quad \gamma \quad {e\left( {s,\beta} \right)}}}} \right)}\quad \frac{\gamma}{\sin \quad \gamma}}}},} \\{{{\Psi_{2}\left( {s,\beta} \right)} = {\int_{0}^{2\quad \pi}{\left( {\nabla_{u{({s,x})}}D_{f}} \right)\left( {{y(s)},{{\cos \quad \gamma \quad \beta} + {\sin \quad \gamma \quad {e\left( {s,\beta} \right)}}}} \right)\quad {\cot (\gamma)}{\gamma}}}},} \\{{\Psi_{3}\left( {s,\beta} \right)} = {\int_{0}^{2\quad \pi}{\left( {\frac{\partial\quad}{\partial\gamma}{D_{f}\left( {{y(s)},{{\cos \quad \gamma \quad \beta} + {\sin \quad \gamma \quad {e\left( {s,\beta} \right)}}}} \right)}} \right)\quad \frac{\gamma}{\sin \quad \gamma}}}}\end{matrix}$

 for all β=β_(j). As before, by itself convolution with a filter is wellknown in the field and can be implemented, for example, as described inU.S. Pat. No. 5,881,123 to Tam, which is incorporated by reference.

Step 43. Using the processed CB data (∇_(u(s,x))D_(f))(y(s),β_(j))evaluate numerically the integralΨ₄(s, β) = ∫₀^(2  π)(∇_(u(s, x))D_(f))(y(s), cos   γ  β + sin   γ  e(s, β))  γ.

Step 44. Store the results of Step 42 and 43 in computer memory.

Step 47. Repeat Steps 41-44 for all lines in the said family of lines.

Step 50. Back-projection

Steps 51-53 are the same as in Embodiment One with s₀ replaced by s.

Step 54. Using equation (9) identify the lines from the said family oflines and points on the said lines that are close to the said projection{circumflex over (x)}. This will give a few values of Ψ_(k)(s,β_(j)),k=1,2,3,4, for β_(j) close to β(s,x).

Step 55. With interpolation compute the quantities Ψ_(k)(s,β(s,x)),k=1,2,3,4, from the said values Ψ_(k)(s,β_(j)), k=1,2,3,4, respectively,for β_(j) close to β(s,x).

Step 56. Multiply the quantities Ψ_(k)(s,β(s,x)), k=1,2,3,4, by thefactors taken from equation (17) and add them up. This will give thequantity${A\left( {s,x} \right)} = {{\left( {\frac{\partial\quad}{\partial s}\frac{1}{{x - {y(s)}}}} \right){\Psi_{1}\left( {s,{\beta \left( {s,x} \right)}} \right)}} + {\frac{{\beta_{s}\left( {s,x} \right)} \cdot {u\left( {s,x} \right)}}{{x - {y(s)}}}{\Psi_{2}\left( {s,{\beta \left( {s,x} \right)}} \right)}} + {\frac{{e_{s}\left( {s,x} \right)} \cdot {u\left( {s,x} \right)}}{{x - {y(s)}}}{\Psi_{4}\left( {s,{\beta \left( {s,x} \right)}} \right)}} + {\frac{{\beta_{s}\left( {s,x} \right)} \cdot {e\left( {s,x} \right)}}{{x - {y(s)}}}{{\Psi_{3}\left( {s,{\beta \left( {s,x} \right)}} \right)}.}}}$

Step 57. Add the said quantity A(s,x) to the image being reconstructedat the point x according to a pre-selected scheme (for example, theTrapezoidal scheme) for approximate evaluation of the integral withrespect to s in equation (17).

Step 58. If the parameter value s corresponding to the CB projection,which is currently in computer memory, is close to a boundary point ofthe parametric interval I_(PI)(x)—either s_(b)(x) or s_(t)(x), thenusing interpolation find$\frac{\Psi_{1}\left( {s^{\prime},{\beta \left( {s^{\prime},x} \right)}} \right)}{{x - {y\left( s^{\prime} \right)}}}.$

 Here s′ is either s_(b)(x) or s_(t)(x). Using the found value updatethe image being reconstructed at the point x according to equation (17).

Step 59. Go to Step 51 and choose a different reconstruction point x.

Step 60.

Step 61. Fix a reconstruction point x. If all the subsequent CBprojections are not needed for reconstructing the image at this point,divide the value of the computed image at x by −2π² and display theimage at x on the computer display 6 of FIG. 1. Repeat this step for allthe reconstruction points.

Step 62. If not all the CB projections have been processed, go to Step10 and load the next CB projection into computer memory. The algorithmconcludes if the remaining CB projections are not needed for imagereconstruction at any of the reconstruction points x or if there are nomore CB projections to process.

Embodiment Three

For the third embodiment, two families of lines can be used incombination with either one CB projection in memory or with several CBprojections in memory. We will first present the main reconstructionformula, and then describe the novel algorithm. Denote $\begin{matrix}{{e_{1}\left( {s,x} \right)}:={\frac{\left\lbrack {{\beta \left( {s,x} \right)} \times {\overset{.}{y}(s)}} \right\rbrack \times {\beta \left( {s,x} \right)}}{{\left\lbrack {{\beta \left( {s,x} \right)} \times {\overset{.}{y}(s)}} \right\rbrack \times {\beta \left( {s,x} \right)}}}.}} & (19)\end{matrix}$

By construction, e₁(s,x) is a unit vector in the plane through y(s) andspanned by β(s,x),{dot over (y)}(s). Moreover, e₁(s,x)⊥β(s,x). Giveny(s),sε(s_(b)(x),s_(t)(x)), {haeck over (s)}(x), finds_(tan)εI_(PI)(x),s_(tan)≠s, such that the plane through x,y(s), andy(s_(tan)) is tangent to C_(PI)(x) at y(s_(tan)). This is equivalent tosolving

[(x−y(s _(tan)))×(x−y(s))]·{dot over (y)}(s _(tan)))=0, s _(tan)≠s.  (20)

Existence and uniqueness of the solution s_(tan)εI_(PI)(x) to (20) isshown in A. Katsevich “Theoretically exact FBP-type inversion algorithmfor Spiral CT”, which is accepted in January 2002 for publication in the“SIAM Journal on Applied Mathematics”. It is also shown there thats_(tan)(s,x) is smooth with respect to s on (s_(b)(x),s_(t)(x)), {haeckover (s)}(x) and is made continuous on [s_(b)(x),s_(t)(x)] by setting

s _(tan)(s,x)=s _(t)(x) if s=s _(b)(x)

s _(tan)(s,x)={haeck over (s)}(x) if s={haeck over (s)}(x),  (21)

s _(tan)(s,x)=s _(b)(x) if s=s _(t)(x).

Once s_(tan)=s_(tan)(s,x) has been found, denote similarly to (19)$\begin{matrix}{{{e_{2}\left( {s,x} \right)}:=\frac{\left\lbrack {{\beta \left( {s,x} \right)} \times \Theta} \right\rbrack \times {\beta \left( {s,x} \right)}}{{\left\lbrack {{\beta \left( {s,x} \right)} \times \Theta} \right\rbrack \times {\beta \left( {s,x} \right)}}}},} & (22)\end{matrix}$

where

Θ=sgn(s−s _(tan)(s,x))β(s _(tan) ,x) if sε(s _(b)(x),s _(t)(x)), {haeckover (s)}(x),

Θ={dot over (y)}(s _(tan)) if sε{s _(b)(x),{haeck over (s)}(x),s_(t)(x)}.  (23)

By construction, e₂(s,x) is a unit vector in the plane through x,y(s),and tangent to C_(PI)(x) at y(s_(tan)). In addition, e₂(s,x)⊥β(s,x). Theinversion formula is now as follows: $\begin{matrix}{\quad {{f = {\frac{1}{2}\left( {{B_{1}f} + {B_{2}f}} \right)}},\quad {where}}} & (24) \\{{\left( {B_{k}f} \right)(x)}:={{- \frac{1}{2\quad \pi^{2}}}{\int_{I_{PI}{(x)}}^{\quad}{\frac{1}{{x - {y(s)}}} \times {{\left. {\int_{0}^{2\quad \pi}{\frac{\partial\quad}{\partial q}{D_{f}\left( {{y(q)},{{\cos \quad \gamma \quad {\beta \left( {s,x} \right)}} + {\sin \quad \gamma \quad {e_{k}\left( {s,x} \right)}}}} \right)}}} \middle| {}_{q = s}\quad {\frac{\gamma}{\sin \quad \gamma}\quad {s}} \right.,{k = {1,2.}}}}}}}} & (25)\end{matrix}$

This formula was proven to be theoretically exact in A. Katsevich“Theoretically exact FBP-type inversion algorithm for Spiral CT”, whichwas accepted in January 2002 for publication in the “SIAM Journal onApplied Mathematics”.

Now we describe an efficient (that is, of the FBP type) implementationof inversion formula (24), (25). Denoting $\begin{matrix}{{{e_{1}\left( {s,\beta} \right)}:=\frac{\left\lbrack {\beta \times {\overset{.}{y}(s)}} \right\rbrack \times \beta}{{\left\lbrack {\beta \times {\overset{.}{y}(s)}} \right\rbrack \times \beta}}},{\beta \in S^{2}},} & (26)\end{matrix}$

rewrite B₁f as follows: $\begin{matrix}{{{\left( {B_{1}f} \right)(x)}:={{- \frac{1}{2\quad \pi^{2}}}{\int_{I_{PI}{(x)}}^{\quad}{\frac{1}{{x - {y(s)}}}{\Psi_{1}\left( {s,{\beta \left( {s,x} \right)}} \right)}{s}}}}},} & (27) \\{{\Psi_{1}\left( {s,\beta} \right)}:=\left. {\int_{0}^{2\quad \pi}{\frac{\partial\quad}{\partial q}{D_{f}\left( {{y(q)},{{\cos \quad \gamma \quad \beta} + {\sin \quad \gamma \quad {e_{1}\left( {s,\beta} \right)}}}} \right)}}} \middle| {}_{q = s}{\frac{1}{\sin \quad \gamma}\quad {{\gamma}.}} \right.} & (28)\end{matrix}$

Let Π(ω),ωεR, where R is the set of real numbers, denote the family ofplanes containing y(s) and parallel to {dot over (y)}(s). Intersectionsof Π(ω) with the detector plane generate a family of lines L(ω) parallelto L₀ (see FIG. 12). Fix any βεΠ(ω). By construction, vectors cos γβ+sinγe₁(s,β),0≦γ<2π, belong to the same plane Π(ω). Recall that forconvenience we think of vectors β,e₁(s,β), and their linear combinationsas if they are attached to y(s). Let θ be a polar angle in Π(ω). Sincee₁(s,β)·β=0,|e₁(s,β)|=1, we can write:

β=(cos θ,sin θ),e ₁(s,β)=(−sin θ,cos θ),β,e ₁(s,β)εΠ(ω).  (29)

Therefore, $\begin{matrix}{{\Psi_{1}\left( {s,\beta} \right)} = {\int_{0}^{2\pi}{\frac{\partial\quad}{\partial q}\quad {D_{f}\left( {{y(q)},\left( {{\cos \left( {\theta + \gamma} \right)},{\sin \left( {\theta + \gamma} \right)}} \right)} \right)}{_{q = s}{{\frac{1}{\sin \quad \gamma}{\gamma}},{\beta \in {{\Pi (\omega)}.}}}}}}} & (30)\end{matrix}$

Equation (30) is of convolution type. Hence, one application of FFT tothe integral in equation (30) gives values of Ψ₁(s,β) for all βεΠ(ω) atonce. Calculation of B₂f can be arranged in a similar way. It followsfrom equation (20) that, apart from the condition s_(tan)εI_(PI)(x),s_(tan) actually depends only on s and β(s,x). Therefore, we can write$\begin{matrix}{{{e_{2}\left( {s,\beta} \right)}:=\frac{\left\lbrack {\beta \times u} \right\rbrack \times \beta}{{\left\lbrack {\beta \times u} \right\rbrack \times \beta}}},{u = {u\left( {s,\beta} \right)}},{\beta \in S^{2}},} & (31) \\{{\Psi_{2}\left( {s,\beta} \right)}:={\int_{0}^{2\pi}{\frac{\partial\quad}{\partial q}\quad {D_{f}\left( {{y(q)},{{\cos \quad {\gamma\beta}} + {\sin \quad \gamma \quad {e_{2}\left( {s,\beta} \right)}}}} \right)}{_{q = s}{{\frac{1}{\sin \quad \gamma}{\gamma}},}}}}} & (32) \\{{\left( {B_{2}f} \right)(x)}:={{- \frac{1}{2\pi^{2}}}{\int_{I_{PI}{(x)}}{\frac{1}{{x - {y(s)}}}{\Psi_{2}\left( {s,{\beta \left( {s,x} \right)}} \right)}{{s}.}}}}} & (33)\end{matrix}$

Fix s_(tan)ε[s−2π+Δ,s+2π−Δ],s_(tan)≠s, and let Π(s_(tan)) denote theplane through y(s),y(s_(tan)), and containing {dot over (y)}(s_(tan)).If s_(tan)=s, Π(s_(tan)) is determined by continuity and coincides withthe plane through y(s) and parallel to {dot over (y)}(s), ÿ(s). Thefamily of lines L(s_(tan)) obtained by intersecting Π(s_(tan)) with thedetector plane is shown in FIG. 13. By construction, given any xεU withβ(s,x)εΠ(s_(tan)) and such that {circumflex over (x)} appears to theleft (right) of the point of tangency s_(tan) if {circumflex over (x)}is above (below) L₀, then s_(tan) used here is precisely the same ass_(tan) provided by equations (20) and (21). The condition that we haveformulated regarding the location of {circumflex over (x)} relative tos_(tan) and L₀ guarantees that s_(tan)εI_(PI)(x). Sincee₂(s,β)·β=0,|e₂(s,β)|=1, we can write:

β=(cos θ,sin θ),e ₂(s,β)=(−sin θ,cos θ),β,e ₂(s,β)εΠ(s _(tan)).  (34)

Therefore, $\begin{matrix}{{{\Psi_{2}\left( {s,\beta} \right)} = \left. {\int_{0}^{2\pi}{\frac{\partial\quad}{\partial q}\quad {D_{f}\left( {{y(q)},\left( {{\cos \quad \left( {\theta + \gamma} \right)},{\sin \left( {\theta + \gamma} \right)}} \right)} \right)}}} \middle| {}_{q = s}{\frac{1}{\sin \quad \gamma}{\gamma}} \right.},{\beta \in {{\Pi \left( s_{\tan} \right)}.}}} & (35)\end{matrix}$

Equation (35) is of convolution type and one application of FFT givesvalues of Ψ₂(s,β) for all βεΠ(s_(tan)) at once. After Ψ₂(s,β) has beencomputed, we use only the portion of it that is located to the left(right) of the point of tangency s_(tan) if L(s_(tan)) is above (below)L₀. In numerical implementation of equations (27), (30), and (33), (35)one can use bilinear interpolation to pass from a rectangular grid ofpoints on the detector to points on the lines L(ω) and L(s_(tan)) andback. As suggested by equations (30) and (35), the points on L(ω) andL(s_(tan)) can be parameterized by polar angle in the correspondingplane.

Equations (27), (30), and (33), (35) demonstrate that the resultingalgorithm is of the FBP type. First, one computes shift-invariantfiltering of a derivative of CB projections using (30) for all requiredω: ω_(min)≦ω≦ω_(max) (cf. FIG. 12), and using equation (35)—for alls_(tan)ε[s−2π+Δ,s+2π−Δ] (cf. FIG. 13). The second step isback-projection according to (27) and (33). Since ∂/∂q in equations (30)and (35) is a local operation, each CB projection is stored in memory assoon as it has been acquired for a short period of time for computingthis derivative at a few nearby points and is never used later. Now wedescribe the algorithm in detail.

Step 10 is the same as Step 10 in Embodiment One.

Step 20. Selecting the two sets of lines.

Step 21. Choose a discrete set of values of the parameter ω inside theinterval ω_(min)≦ω≦ω_(max) (see FIG. 12). This will give a collection ofplanes Π(ω_(j)) containing y(s₀) and parallel to {dot over (y)}(s₀).

Step 22. Intersections of Π(ω_(j)) with the detector plane DP(s₀)generates the first family of lines L(ω_(j)) parallel to L₀ (see FIG.12).

Step 23. Choose a discrete set of values of the parameter s_(tan,j)inside the interval [s₀−2π+Δ,s₀+2π−Δ].

Step 24. For each selected s_(tan)=s_(tan,j) consider the planeΠ(s_(tan)) through y(s),y(s_(tan)), and containing {dot over(y)}(s_(tan)). If s_(tan)=s, Π(s_(tan)) is determined by continuity andcoincides with the plane through y(s₀) and parallel to {dot over(y)}(s₀),ÿ(s₀).

Step 25. The family of lines L(s_(tan,j)) obtained by intersectingΠ(s_(tan,j)), selected in Step 24, with the detector plane DP(s₀) is therequired second family of lines (see FIG. 13).

Step 30. Preparation for filtering

This step is essentially the same as Step 30 of Embodiment One. Theminor differences are as follows. Here this step is used twice. Thefirst time it is applied to the first family of lines L(ω_(j)) and givesthe first processed CB data Ψ₁(s₀,β_(j)). The second time it is appliedto the second family of lines L(s_(tan,j)) and gives the secondprocessed CB data Ψ₂(s₀,β_(j)).

Step 40. Filtering

The filtering step is also essentially the same as Step 40 of EmbodimentOne. The only difference is that here it is used twice. The first timeit is used to convolve the first processed CB data Ψ₁(s₀,β_(j)) withfilter 1/sin γ along lines L(ω_(j)) from the first family of linesgiving the first filtered CB data Φ₁(s₀,β_(j)). The second time it isused to convolve the second processed CB data Ψ₂(s₀,β_(j)) with filter1/sin γ along lines L(s_(tan,j)) from the second family of lines givingthe second filtered CB data Φ₂(s₀,β_(j)).

Step 50. Back-projection

The back-projection step is also essentially the same as Step 50 ofEmbodiment One. The only difference is that here Steps 51-56 are usedtwice. The first time they are used to back-project the first filteredCB data Φ₁(s₀,β_(j)). The second time they are used to back-project thesecond filtered CB data Φ₂(s₀,β_(j)).

Step 57. Following equation (24), add the said contributions from thefirst and second back-projected CB data to the image being reconstructedat the point x according to a pre-selected scheme (for example, theTrapezoidal scheme) for approximate evaluation of the integrals inequation (25).

Step 58 is the same as Step 58 in Embodiment One.

Step 60 is the same as Step 60 in Embodiment One.

Other Embodiments of the invention are possible. For example, one canintegrate by parts in equation (25) (similarly to what was done withequation (10)—see equation (17)), to get an exact FBP-type inversionformula which requires keeping only one CB projection in computermemory. The algorithmic implementation of this alternative embodimentwill be very similar to the algorithmic implementation of EmbodimentTwo.

While the invention has been described, disclosed, illustrated and shownin various terms of certain embodiments or modifications which it haspresumed in practice, the scope of the invention is not intended to be,nor should it be deemed to be, limited thereby and such othermodifications or embodiments as may be suggested by, the teachingsherein are particularly reserved especially as they fall within thebreadth and scope of the claims here appended.

I claim:
 1. A method of reconstructing images from data provided by atleast one detector, comprising the steps of: scanning an object by atleast one cone beam projection and at least one detector; andreconstructing an exact image of the scanned object with a convolutionbased FBP (Filtered Back Projection) algorithm.
 2. The method of claim1, wherein the step of scanning includes the step of: spiral scanningthe object.
 3. The method of claim 1, wherein the scanning step furtherincludes the step of: moving a table supporting the object through arotational scanner.
 4. The method of claim 1, wherein the step ofreconstructing further includes the step of: shift invariant filteringcone beam projections; and back projection updating the image of thescanned object.
 5. The method of claim 1, wherein the step ofreconstructing includes the steps of: storing approximately 2 toapproximately 4 cone beam (CB) projections in memory at a time; andusing one family of lines for the step of reconstructing.
 6. The methodof claim 1, wherein the step of reconstructing includes the steps of:storing 1 cone beam (GB) projection in memory at a time; and using onefamily of lines for the step of reconstructing.
 7. The method of claim1, wherein the step of reconstructing includes the steps of: storingapproximately 2 to approximately 4 cone beam (GB) projections in memoryat a time; and using two families of lines for the step ofreconstructing.
 8. The method of claim 1, wherein the step ofreconstructing includes the steps of: storing 1 cone beam (GB)projection in memory at a time; and using two families of lines for thestep of reconstructing.
 9. A method of computing images derived fromspiral computer tomography with detectors, comprising the steps of: (a)collecting cone beam data from a detector during a scan of an object;(b) identifying lines on a plane Π intersecting the cone beam, whereinthe step (b) of identifying lines includes the steps of: (bi) choose adiscrete set of values of the parameter s₂ inside an interval containings₀, where s₀ and s₂ are values of the parameter along the spiral; (bii)compute the vector u(s₀, s₂) for each selected s₂ according to equations${u\left( {s_{0},s_{2}} \right)} = {\frac{\left( {{y\left( s_{1} \right)} - {y\left( s_{0} \right)}} \right) \times \left( {{y\left( s_{2} \right)} - {y\left( s_{0} \right)}} \right)}{{\left( {{y\left( s_{1} \right)} - {y\left( s_{0} \right)}} \right) \times \left( {{y\left( s_{2} \right)} - {y\left( s_{0} \right)}} \right)}}{{sgn}\left( {s_{2} - s_{0}} \right)}}$${{{if}\quad 0} < {{s_{2} - s_{0}}} < {2\pi}},{{u\left( {s_{0},s_{2}} \right)} = {{\frac{{\overset{.}{y}\left( s_{0} \right)} \times {\overset{¨}{y}\left( s_{0} \right)}}{{{\overset{.}{y}\left( s_{0} \right)} \times {\overset{¨}{y}\left( s_{0} \right)}}}\quad {if}\quad s_{2}} = {s_{0}.}}}$

where, y(s₀), y(s₁), y(s₂) are three points on the spiral relatedaccording to the following equations s ₁=ψ(s ₂ −s ₀)+s ₀ if s ₀ ≦s ₂ <s₀+2π,  (4) s ₁=ψ(s ₀ −s ₂)+s ₂ if s ₀−2π<s ₂ <s ₀.  (5) and ψ is asmooth function with the properties ψ(0)=0; 0<ψ′(t)<1, tε[0,2π]; u(s₀,s₂) is a unit vector perpendicular to the plane containing the pointsy(s₀), y(s₁), y(s₂); {dot over (y)}(s):=dy/ds; ÿ(s):=d²y/ds²; (biii)find a line for each u(s₀, s₂) which is obtained by intersecting theplane through y(s₀) and perpendicular to the said vector u(s₀,s₂) withthe plane Π; and (biv) repeating steps (bi-biii), and forming a familyof lines from a collection of the lines; (c) preprocessing and shiftinvariant filtering said data along said lines; (d) back projecting saidfiltered data to form a precursor of said image; and (e) repeating stepsa, b, c, and d until an image of the object is completed.
 10. The methodof claim 9, wherein the scan includes an x-ray exposure of the object.11. The method of claim 9, further comprising the steps of: storing 1cone beam(CB) projection in memory at a time; and using two families oflines for reconstructing the exact image.
 12. The method of claim 9,further including the step of (b5) preparing the data prior to the step(c) filtering step of: (b5i) fixing a line L(s₂) from the set of linesobtained in step (b); (b5ii) parameterizing points on the lines by polarangle γ in the plane through y(s₀) and L(s₂); (b5iii) choosing adiscrete set of equidistant values γ_(j); (b5iv) find the unit vectorβ_(j) for each γ_(j) which points from y(s₀) towards the point on L(s₂)that corresponds to γ_(j); (b5v) computer derivatives using the GBprojection data D_(f)(y(q),Θ) for a few values of q close to s₀ usingequation (∂/∂q)D_(f)(y(q),Θ|_(q=s) ₀ for all Θ=β_(j); (b5vi) store thecomputed derivatives; and (b5vii) repeat steps (b5i) to b5vi for alllines L(s₂) identified step (b), in order to create processed CB dataΨ(s₀,γ_(j)) corresponding to the x-ray source located at y(s₀).
 13. Themethod of claim 9, wherein the back-projection step (d) includes thesteps of: (di) fix a reconstruction point x, which represents a pointinside the object being scanned where it is required to reconstruct theimage; (dii) if s₀ belongs to I_(PI)(x), then the said filtered CB dataaffects the image at x and one performs Steps (diii) to (dviii), if s₀is not inside interval I_(PI)(x), then the said filtered CB data is notused for image reconstruction at x and go back to step (di) and chooseanother reconstruction point. Here I_(PI)(x) the PI parametric interval;(diii) find the projection {circumflex over (x)} of x onto a detectorplane DP(s₀) and unit vector γ(s₀,x), which points from y(s₀) towards x;(div) identify lines from family of lines and points on the said linesthat are close to the said projection {circumflex over (x)}, usingequation (x−y(s ₀))·u(s₀,s₂)=0,s₂εI_(PI)(x) (dv) using interpolate findvalue of φ(s₀,β(s₀,x)) from φ(s₀,β_(j)) for β_(j) close to β(s₀,x);(dvi) compute contribution from filtered CB data to the image beingreconstructed at the point x by dividing φ(s₀β(s₀,x)) by −2π²|x−y(s₀)|;(dvii) add the contribution to the image being reconstructed at thepoint x according to a pre-selected scheme; and (dviii) go to step (di)and choose a different reconstruction point x.
 14. The method of claim9, further comprising the steps of: storing approximately 2 toapproximately 4 cone beam (CB) projections in memory at a time; andusing one family of lines for the step of reconstructing.
 15. The methodof claim 9, further comprising the steps of: storing 1 cone beam(CB)projection in memory at a time; and using one family of lines forreconstructing the exact image.
 16. The method of claim 9, furthercomprising the steps of: storing approximately 2 to approximately 4 conebeam(CB) projections in memory at a time; and using two families oflines for reconstructing the exact image.
 17. A method of identifyinglines on a plane Π used for reconstructing images based on a spiral scanof an object in a computer tomography system, comprising the steps of:(i) choose a discrete set of values of the parameter s₂ inside aninterval containing s₀, where s₀ and s₂ are values of the parameteralong the spiral scan path; (ii) compute the vector u(s₀, s₂) for eachselected s₂ according to equations${u\left( {s_{0},s_{2}} \right)} = {\frac{\left( {{y\left( s_{1} \right)} - {y\left( s_{0} \right)}} \right) \times \left( {{y\left( s_{2} \right)} - {y\left( s_{0} \right)}} \right)}{{\left( {{y\left( s_{1} \right)} - {y\left( s_{0} \right)}} \right) \times \left( {{y\left( s_{2} \right)} - {y\left( s_{0} \right)}} \right)}}{{sgn}\left( {s_{2} - s_{0}} \right)}}$${{{if}\quad 0} < {{s_{2} - s_{0}}} < {2\pi}},{{u\left( {s_{0},s_{2}} \right)} = {{\frac{{\overset{.}{y}\left( s_{0} \right)} \times {\overset{¨}{y}\left( s_{0} \right)}}{{{\overset{.}{y}\left( s_{0} \right)} \times {\overset{¨}{y}\left( s_{0} \right)}}}\quad {if}\quad s_{2}} = {s_{0}.}}}$

where, y(s₀), y(s₁), y(s₂) are three points on the spiral relatedaccording to the following equations s ₁=ψ(s ₂ −s ₀)+s ₀ if s ₀ ≦s ₂ <s₀+2π, s ₁=ψ(s ₀ −s ₂)+s ₂ if s ₀−2π<s ₂ <s ₀, and ψ is a smooth functionwith the properties ψ(0)=0; 0<ψ′(t)<1, tε[0,2π];u(s₀,s₂) is a unitvector perpendicular to the plane containing the points y(s₀), y(s₁),y(s₂); {dot over (y)}(s):=dy/ds; ÿ:=d²y/ds²; (iii) find a line for eachu(s₀, s₂) which is obtained by intersecting the plane through y(s₀) andperpendicular to the said vector u(s₀,s₂) with the plane Π; and (iv)repeating steps (i-iii), and forming a family of lines from a collectionof the lines so that an image is reconstructed using the lines.